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We introduce tiie Modular Algorithm for Relativistic Treatment of heavy loN Interactions 
(MARTINI), a comprehensive event generator for the hard and penetrating probes in high energy 
nucleus-nucleus collisions. Its main components are a time evolution model for the soft background, 
PYTHIA 8.1 and the McGill- AMY parton evolution scheme including radiative as well as elastic 
processes. This allows us to generate full event configurations in the high pr region that take into 
account thermal QCD and QED effects as well as effects of the evolving medium. We present results 
for the neutral pion nuclear modification factor in Au-|-Au collisions at RHIC as a function of pr 
for different centralities, and also as a function of the angle with respect to the reaction plane for 
non-central collisions. Furthermore, we study the production of high transverse momentum photons 
incorporating a complete set of photon-production channels. 



I. INTRODUCTION 

High transverse momentum jets emerging from the 
central rapidity region in heavy-ion collisions provide im- 
portant information on the produced hot quark-gluon 
plasma (QGP). To extract this information from the ex- 
perimental data, it is important to develop a good theo- 
retical understanding of the interactions of hard partons 
with the medium. Therefore the energy loss of hard par- 
tons traversing the medium has been under extensive the- 
oretical investigation. Gluon bremsstrahlung including 
the Landau-Pomeranchuk-Migdal (LPM) [l[ effect has 
been proposed as the dominant mechanism for energy loss 
and different theoretical formalisms have been developed 
and applied to describe it 0, i, i, i, i, 0, i, i, [13, [IH, 111 ■ 
Also, binary elastic scattering off thermal partons is po- 
tentially important for the energy loss and momentum 
broadening of high-pr partons. I n (1311 . elastic energy 
loss was combined with the AMY |lOl. Illl. [T^ radiative 
energy loss within the McGill- AMY evolution formalism, 
and in |lj| the description of the coUisional processes was 



further improved. 

The main goal of these and other calculations is to 
create a quantitative basis for the "tomography" of the 
hot and dense nuclear medium created in heavy-ion col- 
lisions. Among the "tomographic variables" is the nu- 
clear modification factor Raa, which is defined as the 
ratio of the hadron yield in A-l-A collisions to that in 
binary-scaled p-|-p interactions. A variety of computa- 
tions that differ significantly in the applied energy loss 
mechanism can reproduce the measured Raa in central 
Au-hAu colhsions at -/s = 200 AGeV at the Relativistic 
Heavy-Ion Collider (RHIC), given the present experimen- 
tal errors (see [isj for a systematic comparison of three 
different formalisms). For this reason the tomographic 
usefulness of Raa in central collisions has become ques- 
tionable More differential observables, such as Raa 
as a function of both pT and azimuth in non-central col- 
lisions [13, [H, [13, HO] J can help to improve this situation. 

Other important observables in heavy-ion collisions are 
electromagnetic probes, such as photons and dileptons. 
Because of the small electromagnetic coupling, once pro- 



duced, they usually escape the medium without further 
interaction and thus provide undistorted information on 
the early stages of a heavy-ion collision. Photon pro- 
duction from nuclear collisions at RHIC has been cal- 
culated in [21I, [2^, [2^ using 1-1-1 dimensional Bjorken 
evolution, 2-1-1 dimensional, or 3-1-1 dimensional relativis- 
tic hydrodynamic evolution for the background, respec- 
tively. Photons in nuclear collisions come from a variety 
of sources, namely direct photons, fragmentation pho- 
tons, jet-plasma photons and thermal photons. Direct 
photons are predominantly produced from hard collisions 
in the early stage of the heavy-ion collision via annihila- 
tion and Compton scattering processes. Fragmentation 
photons are produced from the surviving high energy jets 
after their passing through the thermal medium. Ther- 
mal photons have a negligible contribution at high pT 
and thus are not relevant in this range. On the other 
hand, jet-plasma photons from photon radiation q — > (77 
(bremsstrahlung photons) and jet-photon conversion via 
Compton and annihilation processes have been shown to 
be very important for the understanding of experimental 
data for photon production in Au-|-Au collisions at RHIC 
[Il|,[13,[2l,|2l. 

For a best possible comparison of the theoretical de- 
scription with experimental data, we incorporate the 
McGill- AMY formalism [l^, [HI, [3, [H, [pi for radiative 
energy loss as well as elastic processes [ij into a new 
event generator, dubbed Modular Algorithm for Rela- 
tivistic Treatment of heavy loN Interactions (MARTINI) . 
Its main ingredients are PYTHIA 8.1 [2i,[23 to gener- 
ate the hard partons from the individual nucleon-nucleon 
collisions and to take care of the final fragmentation into 
hadrons, the evolution of the background medium (e.g. 
from hydrodynamic models), and the McGill- AMY par- 
ton evolution formalism. This way it is possible to study 
hard observables in heavy-ion collisions theoretically on 
an event-by-event basis, keeping information on correla- 
tions. 

Several Monte Carlo simulations for heavy-ion colli- 
sions have been or are being developed [H, [2^, [13, [Ml, 
[321 . [33I . [34I . [35I . [36j . The implementation of medium ef- 
fects and use of approximations varies significantly be- 
tween the different models. MARTINI is the first Monte 
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Carlo simulation to include AMY bremsstrahlung com- 
bined with elastic processes and in-medium photon pro- 
duction. It is very flexible due to its ability to incor- 
porate any soft background evolution that provides in- 
formation on the temperature and flow velocities of the 
medium. Currently, hydrodynamic evolution calcula- 
tions from four different groups [13, [H, [H, SO, \M, El 
have been implemented. 

In this work we present first results on hard observables 
obtained with MARTINI. These include the azimuth av- 
eraged Raa as a function of px, Raa as a function of 
the azimuth and pr, as well as photon yields and photon 
Raa- 

This paper is organized as follows. We review the 
McGill-AMY evolution formalism and present transition 
rates for radiative and elastic processes in Section|TTl The 
Monte Carlo simulation is introduced in Section HITl re- 
sults are presented in Section IIVI and we conclude and 
present an outlook in Section |Vl 



They are given by [2l|, |43 



II. THE MCGILL-AMY EVOLUTION 
FORMALISM 



At the core of MARTINI lies the McGill-AMY formal- 
ism for jet evolution in a dynamical thermal medium. 
Here, we briefly review this formalism and discuss its im- 
plementation into MARTINI in the following section. 

In previous works [l^, [H, the evolution of the jet 
momentum distribution in a hydrodynamic background 
was computed within the McGill-AMY scheme. This evo- 
lution is governed by a set of coupled Fokker-Planck type 
rate equations of the form 
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where g is the strong coupling constant (as = (7^/(47r)), 
A^f is the number of flavors, Cs is the quadratic Casimir 
relevant for the process, and x = k/p is the momentum 
fraction of the radiated gluon (or the quark, for the case 
g qq). h = (k x p) x e|| determines how non-collinear 
the final state is, where e|| is the unit vector in a direction 
nearly coUincar with k, p and p -f- k that can be fixed 
by convention (see Parametrically h is of order 

gT^ and hence small compared to p • k. Therefore h 
can be taken as a two-dimensional vector in transverse 
space. F(h, p , k ) is the solution of the following integral 
equation |2ll. l43j: 



2h iSE{h,p,k)F{h 
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x{(a - CA/2)[F(h) - F(h-fcqi)] 

+ (CA/2)[F(h)-F(h+pqi)] 

+ (CA/2)[F(h)-F(h-(p-fc)qx)]}. (3) 

Here 5E{h,p, k) is the energy difference between the final 
and the initial states: 
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and are the medium induced thermal masses. C(qx) 
is the differential rate to exchange transverse (to the || 
direction) momentum qj^. In a hot thermal medium, its 
value at leading order in ag is (4^ 



where dT{p,k)/dk is the transition rate for processes 
where partons of energy p lose energy k. The fc < part 
of the integration incorporates processes which increase 
a parton's energy. 

In the AMY finite temperature field theory approach, 
the energy loss of hard partons is considered in a medium 
at an asymptotically high temperature so that the QCD 
coupling is weak. In this regime, an analytical treatment 
is possible, and a hierarchy of parametrically separated 
scales T > gT > g^T emerges. This allows radiative 
transition rates to be calculated by means of integral 
equations , which correctly reproduce both the Bethe- 
Heitler and the LPM results in their respective limits. 
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■{2N,+Ni). (5) 



For the case of g — > qq, {Cg — Ca/2) should appear as 
the prefactor on the term containing F(h — pqx) rather 
than F(h — fcqx). 

Transition rates for elastic processes can also be com- 
puted in pcrturbative QCD. They are given by [ij] 

du ' V (2^)3 J (2^)3 mpp'qq' 

X 5{p — p' — uj)5{q' — q — uj) 
x\M\y{q,T){l±fiq\T)), (6) 
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where p = IpI and p' = |p'| are the absolute val- 
ues of the three-momenta of the incoming and outgo- 
ing hard parton, respectively, and g = |q| and q' ~ |q'| 
are those of the incoming and outgoing thermal parton. 
to = p—p' = q' — q is the transferred energy. The distribu- 
tion functions / are either Fermi-Dirac or Bose-Einstein 
distributions depending on the nature of the thermal par- 
ton involved. The + or — sign appears accordingly, with 
— for Pauli blocking and + for Bose enhancement, and 
dq describes the degeneracy of the thermal parton. In- 
tegrations over q and q' in Eq. ([6]) can be rewritten in 
terms of the transferred momentum k. Then, employing 
an appropriate separation of the regimes of soft and hard 
momentum exchange, and doing a hard thermal loop re- 
summation in the soft regime to cure the infrared diver- 
gence allows for a numerical evaluation of Eq. It is 
even possible to extend the expression for the soft regime 
to all momenta [3, \^ , which leads to a good approxi- 
mation of the total result and avoids the introduction of 
an intermediate scale. The transition rate as a function 
of both the transferred energy and momentum dT/dujdk 
is obtained from Eq. ([6]) after rewriting and omitting 
the integration over the transferred momentum k. Us- 
ing dT/dwdk will allow us to sample both the transferred 
energy and the transferred three-momentum in an elastic 
process. 

Both radiative and elastic transition rates were com- 
puted and tabulated as functions of the parton energy p, 
as well as the radiated energy k or transferred energy uj 
and momentum k, respectively, and are sampled within 
MARTINI. Note that since the radiated partons are on- 
shell iu = k, while the transferred parton in an elastic 
collision can be far off-shell, so that uj and k are treated 
separately. As opposed to the radiative transition rates, 
the dependence on the coupling of the elastic rates is 
non-trivial and the rates are also tabulated as functions 
of the coupling as- 

Apart from the processes described above we include 
gluon-quark and quark-gluon conversion due to Compton 
and annihilation processes, as well as the QED processes 
of photon radiation q q'y and jet-photon conversion. 
The transition rate for the photon radiation process is 
calculated analogously to the gluon radiation [l^, [l^ . 
and in the limit E ^ T the transition rate for quark- 
photon conversion is given by 



that in Eq. ([7]) times a color factor Cp = 4/3, and 
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where ruq is the thermal quark mass, = g^T'^ /Q, e/ 
is the charge of a quark with flavor /, is the fine 
structure constant, and the delta function indicates that 
we neglect the energy loss of the quark during the process 
of conversion. The conversion rate dr^^/dtj is equal to 



jpconv 

du 



Nf 



m - 1 doj 
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which follows from interchanging the initial with the final 
jet. 



III. MONTE CARLO SIMULATION 

MARTINI solves the rate equations ([T|) using Monte 
Carlo methods. Instead of evolving a probability distri- 
bution P{p), every event and every hard parton is treated 
individually. This way we obtain information on the full 
microscopic event configuration in the high momentum 
regime, including correlations, which allows for a very de- 
tailed analysis and offers a direct interface between the- 
ory and experiment. The average over a large number of 
events will correspond to the solution found by solving 
([l]) for the probability distribution. 

We follow the evolution of one event to describe the 
basic functionality of the simulation. First, the number 
of individual nucleon-nuclcon collisions that produce par- 
tons with a certain minimal transverse momentum p™'" 
is determined from the total inelastic cross-section, pro- 
vided by PYTHIA. The initial transverse positions rj_ of 
these collisions are determined by the initial jet density 
distribution 'PAB{b,r±), which for A-l-B collisions with 
impact parameter b is given by 



TA(r±+h/2)TB{r^-h/2) 

TAB{b) 



(9) 



Here we use a Woods-Saxon form for the nuclear den- 
sity function, p{rj_,z) = po/[l + exp(^^-^)], to evaluate 
the nuclear thickness function TA(rj^) = J dzpAiJ^±, z) 
and the overlap function of two nuclei TAsib) = 
/ d'^r±TA{ri_)TBirj_ + b). The values of the parame- 
ters R = 6.38 fm and d = 0.535 fm are taken from 
The initial parton distribution functions can be selected 
with the help of the Lcs Houches Accord PDF Interface 
(LHAPDF) |47|. Wc include nuclear effects on the parton 
distribution functions using the EKS98 [Hi or EPS08 Hi) 
parametrization, by user choice. In practise we sample 
the initial positions of nucleons in nucleus A and B, su- 
perimpose the transverse areas depending on the impact 
parameter and then divide the overlap region into circles 
with area cTinoi- In three dimensions these are tubes and 
we determine how many jet events with given p™'" occur 
within such a tube using the number of nucleons from 
A and B in the tube, and the probability for a jet event 
jet (Pt™) / '^inci for each of their combinations. 
The soft medium is described by hydrodynamics or 
other models, which provide information on the sys- 
tem's local temperature and flow velocity /3. Before this 
medium has formed, i.e., before the hydrodynamic evo- 
lution begins (t < tq), the partons shower as in vac- 
uum. Currently wc include the complete vacuum shower. 
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because there is no apparent reason why the vacuum 
sphttings should end immediately once the medium has 
formed. Since most of the vacuum shower occurs be- 
fore the medium has formed, this is a reasonable approx- 
imation. We have also tested the other extreme case 
where we stop the vacuum evolution at the virtuality 
scale Qmin = yJvT /tq, determined by the time tq at 
which the medium evolution begins. This modifies the 
parton distribution such that the strong coupling con- 
stant has to be chosen approximately 10% larger to de- 
scribe the pion Raa- At this stage, we do not include 
final state radiation (FSR) of the partons that have left 
the medium - all vacuum showers take place before the 
medium evolution and there is no interference between 
the medium and vacuum showers. The improvement of 
this point is the subject of future work. Once t > rg for a 
certain parton. its in-medium evolution begins. The par- 
tons move through the background according to their ve- 
locity. To compute probabilities for any one of the above 
described in-medium processes, we perform a Lorentz- 
boost into the rest-frame of the fluid cell at the parton's 
position and determine the transition rates according to 
the local temperature and the parton's energy in this lo- 
cal rest-frame. The probability for a parton to undergo 
any process during a time step of length At is given by 
Atr(p,T)totai, where r(p,T)totai = / dkdTtotai/dk for 
the radiative processes, and the integral over both lo and 
k for the elastic processes. dTtotsd/dk is the sum over all 
possible transition rates, which include those for photon 
production. 

In case that some process occurs, we decide which one 
it is according to the relative weights of the different pro- 
cesses at the given temperature and parton energy. We 
sample the radiated or transferred energy from the tran- 
sition rate of the occurring process using the rejection 
method. In case of an clastic process, we also sample 
the transferred transverse momentum, while for radiative 
processes we assume coUinear emission for now, which is a 
good approximation in the weak coupling limit since the 
emission angle is of order g [l^ . After energy and mo- 
mentum are transferred, we boost back into the labora- 
tory frame, where the parton continues to move along its 
trajectory. Radiated partons are also further evolved if 
their momentum is above a certain threshold Pmin, which 
can be set (typically we choose p,nin — 3GeV). This 
leads to a growing and broadening in-medium shower. 
The overall evolution of a parton stops once its energy in 
a fluid cell's rest frame falls below the limit of AT, where 
T is the local temperature. As both radiative and elas- 
tic transition rates were computed in the limit E ^ T, 
it would not make sense to attempt to learn anything 
about the evolution of partons with £' ~ T within this 
formalism. For partons that stay above that threshold, 
the evolution ends once they enter the hadronic phase of 
the background medium. In the mixed phase, processes 
occur only for the QGP fraction. 

When all partons have left the QGP phase, hadroniza- 
tion is performed by PYTHIA, to which the complete 
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FIG. 1: (Color online) Neutral pion spectrum in p-l-p and 
— 10% central Au+Au collisions at RHIC energy compared 
to data by PHENIX [li,[53|. 



information on all final partons is passed. Because 
PYTHIA uses the Lund string fragmentation model 
[51I . [5^ , it is essential to keep track of all color strings 
during the in-medium evolution. This includes gener- 
ating new strings when a gluon is emitted or during a 
conversion process. In the latter case, new string ends 
have to be attached to thermal partons, whose momenta 
arc sampled from Fermi or Bose distributions. Being pro- 
vided with such consistent information on the color string 
structure, PYTHIA is able to perform the fragmentation. 
This concludes the evolution of one event. 

The concept of MARTINI is modular, such that we can 
turn on and off different processes independently, and use 
different hydrodynamic or other data inputs. In principle 
it is also possible to extend MARTINI to use different 
energy loss formalisms. The parameters are set using the 
same interface as PYTHIA and all options for PYTHIA 
can still be modified by the user. 



IV. RESULTS 

A. Pion production and nuclear modification factor 

tv°Raa{pt) 

We begin by showing the MARTINI results for the 
spectrum of neutral pions in p+p collisions at RHIC 
energy {yG = 200 GeV) as well as in central (0-10%) 
Au+Au collisions compared to data by PHENIX [1^, 
[50| in Fig. [1] The calculations were performed using 
CTEQ5L parton distribution functions |53| including nu- 
clear shadowing effects using the EKS98 parametrization 
jist . In the shown results for pion production, isospin 
effects were ignored. We checked that they have no big 
effect (less than 5%), however, we included them in calcu- 
lations of photon production where they become impor- 
tant (15-20%). Au+Au calculations take into account 
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FIG. 2: (Color online) The neutral pion nuclear modification 
factor for mid-central (upper panel) and central (lower panel) 
collisions at RHIC with ^ = 200 GeV. MARTINI resuhs 
(6 — 2.4 fm and b = 7.5 fm) using 3+1 dimensional hydro 
evolution [4l[ and only radiative processes (dashed lines) and 
both radiative and elastic processes (solid lines) with different 
as compared to PHENIX data from [si ]. 



FIG. 3: (Color online) The neutral pion nuclear modification 
factor for mid-central (upper panel) and central (lower panel) 
collisions at RHIC with = 200 GeV. MARTINI resuhs 
(6 — 2.4 fm and b = 7.5 fm) using 2-1-1 dimensional hydro 
evolution [3^ . |39| . |40| | and both radiative and elastic processes 
compared to PHENIX data from 



both radiative and elastic processes in the medium de- 
scribed by either the 2+1 dimensional hydrodynamics of 
3S, [3^, [iOl or the 3+1 dimensional hydrodynamics of 
4ll |. using a coupling constant as = 0.33 or as = 0.3, 
respectively, as was adjusted to describe the experimen- 
tal measurement of the neutral pion nuclear modification 
factor Raa in most central collisions (see below). For 
both hydro evolutions tq = 0.6fm/c. PYTHIA param- 
eters were tuned to fit the p+p data, however there is 
still freedom and possibly an even better set of param- 
eters can be found. Using the same PYTHIA settings, 
we find very good agreement for both p+p and Au+Au 
spectra. Note that in order to achieve higher statistics 
we used certain cutoffs for the minimal px produced in a 
nucleon-nucleon collision in PYTHIA for the Au+Au cal- 
culation. The shown result is a combination of runs with 
minimal p??'" = 1.5 GeV for px < 5.5 GeV, p^"" = 2 GeV 
for 5.5 GeV < pr < 6.5 GeV, and p^'"" = 4 GeV above 
that. Note that the shown results are averages over typ- 
ically ~ 5 • 10® simulated events. 

In Fig. O we present the results for the neutral pion 
nuclear modification factor, defined by 
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1 dNAA{b)/cPpTdy 
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(10) 



in Au+Au collisions at RHIC measured at mid-rapidity 
in two different centrality classes (0-10%) and (20-30%), 
employing the corresponding average impact parameters, 
2.4 fm and 7.5 fm. We compare with the experimental 



measurements by PHENIX [5J]. All parameters are the 
same as in the calculation shown in Fig . [U i.e., hydro- 
dynamic background evolution from |4l| and — 0.3. 
The same value of as is used in all following calculations. 
(There are no additional parameters for the later calcula- 
tion of the azimuthal dependence of Raa or the high-pr 
photon production.) We find very good agreement with 
the data for both centrality classes. 

The value for for which the data is described if we 
only include radiative processes is as — 0.39, higher than 
for the calculations in [ij], where it was 0.33. Differences 
in the treatment of vacuum showers and the fragmenta- 
tion scheme are responsible for this difference. We have 
observed that when varying parameters in PYTHIA, such 
as the factorization or renormalization scale, or param- 
eters in the fragmentation function, we need slighty dif- 
ferent as to describe the data. This freedom leads to an 
approximately 10 - 20% uncertainty in ag. The value 
of as is again closer to that used in [l^ (a^ = 0.27) 
when including elastic processes, because the elastic en- 
ergy loss in our formulation is slightly larger than that 
implemented in [l3| - see [l3l for details on this. Like 
previous works (e.g. [l^) we find that elastic processes 
are important for the computation of Raa- 

Importantly, our approach has the potential to reveal 
the effect of different soft background models on hard 
observables. Fig. [3] shows the result for Raa using the 
2+1 dimensional hydrodynamic evolution from [38l . [39l . 
|40| . As initial temperatures are typically lower than in 
the 3+1 dimensional hydro evolution, we need to increase 
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tts to 0.33 to describe the neutral pion Raa- Further 
effects are demonstrated in the following section. 



B. Angular dependence of the nuclear modification 
factor n^' RaaIpt , ^4>) 

In Figs. [2] and [3] the neutral pion Raa at midrapid- 
ity is averaged over the azimuthal angle with respect to 
the reaction plane (j)- To learn more about the produced 
medium, e.g. to disentangle the effects of the density 
of the medium and the pathlength traversed on the en- 
ergy loss, Raa has been studied at midrapidity in non- 
central collisions not only as a function of pt but also 
as a function of the azimuth (jj in [3, [l^, [lo] ■ We per- 
form the same analysis within MARTINI and separate 
the azimuth in 6 bins of 15° each, reaching from — 15° 
(in-plane) to 75 — 90° (out-of-plane) , as shown in Fig. ID 
Raa is then determined in each bin separately. We em- 
ploy both 2-1-1 and 3-1-1 dimensional hydrodynamic back- 
grounds and ttg = 0.33 and as = 0.3, respectively. Fig. [5] 
shows the most extreme cases (0 — 15° and 75 — 90°) at 
b = 7.5 fm compared to recent mid-central (20-30%) data 
by PHENIX [ifl]. 



out of plane 




FIG. 4: (Color online) Binning in azimuth Acji for the calcu- 
lation of Raa{pt, A<f)) in non-central collisions. 



Particularly at lower px, we find a less pronounced 
difference between the in-plane and out-of-plane result 
than the data, as do previous theoretical calculations 
[T9I . [sst . Furthermore, one can see that the result is 
closer to the experimental data for the 3-f 1 dimensional 
hydrodynamic background. This is due to a larger ini- 
tial eccentricity in this case. Generally, we do not expect 
to describe the low pr part correctly, because we do not 
incorporate soft physics such as flow effects or recombi- 
nation. The too small difference for pt 4 GeV found 
in the simulation could in part stem from a too small 
eccentricity of the initial state when using Glauber ini- 
tial conditions in the hydrodynamic calculation, or a too 
large surface bias in the theoretical description. Use of 
Color-Glass-Condensate initial conditions in 3+1 dimen- 
sional hydrodynamics may help to improve the desciption 
of the data [M, [H^I . 
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FIG. 5: (Color online) Nuclear modification factor as a func- 
tion of Pt and A<^ for mid-central collisions at -^s = 200 GeV. 
MARTINI results (h = 7.5 fm, upper panel 2-l-lD hydro, 
Qfs = 0.33; lower panel 3-l-lD hydro, Os = 0.3) with both 
radiative and elastic processes for most in-plane and most 
out-of-plane angular bins compared to PHENIX data from 
[20I . The gray band indicates the experimental error in Raa ■ 



C. Photon production 

Photons in the high p^ region produced in nuclear 
collisions are dominantly direct photons, fragmentation 
photons, and jet-plasma photons. Direct photons are 
included in PYTHIA. Apart from leading order direct 
photons, PYTHIA produces additional photons emitted 
during the vacuum showers. Parts of these overlap with 
effects that are found in NLO calculations [6l|, [63| , but 
there is no simple theoretical way to identify the amount 
of overlap between the two. Also fragmentation pho- 
tons are part of the photons produced in the showers 
in PYTHIA. We computed photon production within 
PYTHIA with and without photons from showers and 
found that the shower contribution leads to an effective 
if-factor of approximately 1.8 in the regarded pT range. 

For heavy-ion reactions MARTINI adds the very rele- 
vant jet-medium photons from photon radiation 9 — > (77 
(bremsstrahlung photons) and jet-photon conversion via 
Compton and annihilation processes. As above, we in- 
clude the full vacuum shower in the shown calculation - 
most of the shower photons will be emitted before the 
medium has formed and the parton has realized that it 
has formed. To improve on this approximation and in- 
clude interference between vacuum and in-medium radi- 
ation is the subject of future work. 

We present results for photon production in p-l-p and 
Au+Au collisions at -y/s = 200 GeV compared to data by 
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FIG. 6: (Color online) Photon yield in p+p and Au+Au col- 
lisions at RHIC energy ^ = 200 GeV. MARTINI results 
{b = 2.4 fm, 3+lD hydro) compared to data from [SSl. [SOl, [60j . 



FIG. 7: (Color online) Photon nuclear modification factor 
for central collisions at Js = 200 GeV. MARTINI resuhs 
compared to data from [6y] . See main text for details. 



PHENIX [5i,[5i,[63 in Fig. El 

Another observable that provides information on the 
effect of the nuclear medium on photon production is the 
photon nuclear modification factor 

1 dNl^{h)/<PpTdy 
iVcoii(6) dN^p/d^pTdy • ^ ' 

Fig. [7] shows R\j^ as a function of pt for most central 
Au+Au collisions (6 = 2.4 fm) at RHIC compared with 
- 10% central PHENIX data. We find that in the ap- 
proximation of including all photons from the vacuum 
shower leads to good agreement with the data within the 
error bars. We checked that not including any photons 
from the final state vacuum shower before tq, which cor- 
respond to a pre-equilibium contribution, reduces Raa 
by approximately 20%. 



V. CONCLUSIONS AND OUTLOOK 

We presented first results obtained with the newly de- 
veloped Modular Algorithm for Relativistic Treatment 
of heavy loN Interactions (MARTINI). This hybrid ap- 
proach describes the soft background medium using hy- 
drodynamics or other medium models and simulates the 
hard event microscopically, using PYTHIA 8.1 to gener- 
ate individual hard nuclcon-nuclcon collisions. Hard par- 
tons are evolved through the medium using the McGill- 
AMY evolution scheme including radiative and elas- 
tic processes. Fragmentation is performed employing 
PYTHIA 8.1 which uses the Lund string fragmentation 
model. Apart from parameters in PYTHIA which were 
fixed by matching the neutral pion and photon spectra 
in p+p coUisons to experimental data, as is the only free 
parameter. Employing the 3+1 dimensional hydrody- 
namic evolution from it was set to ag = 0.3 to 



match the neutral pion Raa measurement for central col- 
lisions. Using the same value for all other calculations 
(there was no additional freedom in any of the calcula- 
tions), we were able to describe the neutral pion Raa in 
mid-central collisions as well as photon yields and photon 
Raa- For photon production, we added photons from jet- 
medium interactions to the prompt and shower photons 
that PYTHIA produces. 

We have demonstrated that MARTINI is able to re- 
veal the effect of using different evolution models for the 
soft background on results for hard observables. For ex- 
ample, when using a 2+1 dimensional hydrodynamical 
evolution model [H, [3^, , a slightly different value for 
as (tts = 0.33) had to be employed to describe the ex- 
perimental data. The calculated azimuthal dependence 
of the neutral pion Raa in mid-central collisions turned 
out to be weaker than that found by PHENIX, for both 
shown hydrodynamic models, but more so for the 2+1 
dimensional one. This is mainly due to a larger initial 
eccentricity in the 3+1 dimensional model. The result for 
the 3+1 dimensional hydrodynamic background is along 
the line with previous calculations using different energy 
loss schemes |l9l . [Hsf . The weak azimuthal dependence, 
particularly for ^ 6 GeV, can be due to the lack of 
soft physics effects such as flow and recombination, as 
well as a too small eccentricity of the initial state when 
using Glauber initial conditions for the hydro evolution 
calculation. 

Having shown that MARTINI can reproduce one-body 
observables in good agreement with the data, the next 
step will be to explore its full potential by studying many- 
body observables and correlations. Another future task 
is the implementation of heavy quark evolution. 
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